home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / tests / scilab.tst < prev    next >
Text File  |  1999-09-16  |  13KB  |  437 lines
  1. pi=%pi
  2. i=%i
  3. e=%e
  4. //       tests
  5. //
  6. 1
  7. a=1
  8. b=[1 2 3]
  9. c=[1 2 3;4 5 6]
  10. d=[1 2 3;4 5 6]'
  11. d=[1 2 3 4;4 5 6 7;8 9 10 11;12 13 14 15]
  12. d=[1 2;3 4]
  13. e1=[[1 2] [3 4];[5 6 7 8];[9;10;11;12]']
  14. f([1 3 5],[1 2 3])=[-1 -2 -3;-4 -5 -6;-7 -8 -9]
  15. g=[i,2,3;1 i 3;1 2 i]
  16. h=[i 1 2 i 3]
  17. o=[1 i i 2 3]'
  18. b(2)=3
  19. e1(1,3)=1
  20. p='apcdefghijklmnopqrstuvwxyz0123456789'
  21. // test de stackg
  22. a,b,c,d,e1,f,g,h,o,p
  23. b(2)
  24. e1(1,3)
  25. e1([1 2],[3 4])
  26. e1(1:2,:)
  27. e1(:,4)
  28. e1(:,:)
  29. g(1,1)
  30. g(:,1)
  31. g(1:2,:)
  32. g(:,:)
  33. o'
  34. //test de stack2
  35. //
  36. // additions...
  37. [1 2 3;4 5 6]-[1 2 3;4 5 6]
  38. [1 2 3;4 5 6]+2*[1 2 3;4 5 6]-[1 2 3;4 5 6]*3
  39. [1 i 3;4 5 i]+2*[1 i 3;4 5 i]-[1 i 3;4 5 i]*3
  40. [2*i -4*i 2*i]+2*[i -2*i i]+i*[2 -4 2]-[i -2*i i]*2 -2*[2 -4 2]*i
  41. //
  42. 2\[4 8 16]/2-[1 2 4]
  43. i\[2*i 2 i*4]/i-[-i*2 -2 -4*i]
  44. // element wise operations
  45. [1 2;3 4;5 6].*[1 2;3 4;5 6]-[1 4;9 16;25 36]
  46. [1 2;3 4;5 6*i].*[1 2;3 4;5 6]-[1 4;9 16;25 36*i]
  47. [1 2;3 4;5 6].*[1 2;3 4;5 6*i]-[1 4;9 16;25 36*i]
  48. [1 2;3 4;5 6*i].*[1 2;3 4;5 6*i]-[1 4;9 16;25 -36]
  49. //
  50. [2 9 8;3 10 15]./[2 3 4;3 5 5]-[1 3 2;1 2 3]
  51. [2 9*i 8;3 10 15]./[2 3 4;3 5 5]-[1 3*i 2;1 2 3]
  52. [2 9 8;3 10 15]./[2*i 3 4;3 5 5]-[-i 3 2;1 2 3]
  53. [2*i 9 8;3 10 15]./[2*i 3 4;3 5 5]-[1 3 2;1 2 3]
  54. //
  55. [2 3 4;3 5 5].\[2 9 8;3 10 15]-[1 3 2;1 2 3]
  56. [2 i 4;3 5 5].\[2 9 8;3 10 15]-[1 -9*i 2 ;1 2 3]
  57. [2 3 4;3 5 5].\[i 9 8;3 10 15]-[.5*i 3 2;1 2 3]
  58. [i 3 4;3 5 5].\[i 9 8;3 10 15]-[1 3 2;1 2 3]
  59. //multiplication
  60. [1 2 3;4 5 6]*[3;2;1]-[10;28]
  61. [1 i 3;4 5 6]*[3;2;1]-[6+2*i;28]
  62. [1 2 3;4 5 6]*[i;2;1]-[i+7;4*i+16]
  63. [1 i 3;2*i,-i,1]*[i;-i;i]-[4*i+1;-3+i]
  64. //eye
  65. eye(4,4)
  66. //a+-b*eye a*eye+-b a+-eye*b eye*a+-b
  67. -2*eye+[1 2;3 4]+eye*2-[1 2;3 4]
  68. -2*i*eye+[1 2;3 4]+eye*2*i-[1 2;3 4]
  69. // :
  70. 1:10
  71. 1:.1:2
  72. // for
  73. for k=1:3,for l=1:2,a(k,l)=k+l;end;end;a
  74. diag([1 2 3])-[1 0 0;0 2 0;0 0 3]
  75. diag([1 i 2])-[1 0 0;0 i 0;0 0 2]
  76. a=[1 2 3 4;5 6 7 8];
  77. c=a;
  78. c(1,1)=i;
  79. diag(a)-[1;6]
  80. diag(a,1)-[2;7]
  81. diag(a,-1)-[5]
  82. diag(a,4)
  83. diag(c)-[i;6]
  84. diag(c,1)-[2;7]
  85. diag(c,-1)-[5]
  86. diag(c,4)
  87. //
  88. eye(a)
  89. eye(c)
  90. eye(3,3)
  91. eye(2,3)
  92. //
  93. ones(a)
  94. ones(c)
  95. ones(3,3)
  96. ones(3,2)
  97. //
  98. rand(a)
  99. rand(c)
  100. rand(3,3)
  101. rand(3,2)
  102. rand
  103. rand('uniform')
  104. rand('normal')
  105. rand('seed',5)
  106. rand('seed',0)
  107. //
  108. abs(c)
  109. abs(a)
  110. //
  111. real(c)
  112. imag(c)
  113. real(a)
  114. imag(a)
  115. //
  116. round(rand(3,3))
  117. //
  118. conj(a)
  119. conj(c)
  120. //
  121. size(a)
  122. [m,n]=size(a)
  123. //
  124. triu(a)
  125. tril(a)
  126. triu(a,1)
  127. triu(a,-1)
  128. tril(a,1)
  129. tril(a,-1)
  130. triu(c)
  131. tril(c)
  132. triu(c,1)
  133. triu(c,10)
  134. triu(c,-1)
  135. triu(c,-10)
  136. tril(c,1)
  137. tril(c,10)
  138. tril(c,-1)
  139. tril(c,-10)
  140. //test de matlu
  141. a=rand(4,4);b=rand(5,4);ac=a+i*rand(4,4);bc=b+i*rand(5,4);
  142. //
  143. if abs((1/a)*a-eye)> 10*%eps then pause,end
  144. if abs((i/a)*a-i*eye)> 10*%eps then pause,end
  145. if abs((1/ac)*ac-eye)> 10*%eps then pause,end
  146. if abs((i/ac)*ac-i*eye)> 10*%eps then pause,end
  147. if abs(a*(a\1)-eye)> 10*%eps then pause,end
  148. if abs(a*(a\i)-i*eye)> 10*%eps then pause,end
  149. if abs(ac*(ac\1)-eye)> 10*%eps then pause,end
  150. if abs(ac*(ac\i)-eye*i)> 10*%eps then pause,end
  151. //
  152. if abs(inv(a)*a-eye)> 10*%eps then pause,end
  153. if abs(inv(ac)*ac-eye)> 10*%eps then pause,end
  154. //
  155. if abs((b/a)*a-b)> 10*%eps then pause,end
  156. if abs((b/ac)*ac-b)> 10*%eps then pause,end
  157. if abs((bc/a)*a-bc)> 10*%eps then pause,end
  158. if abs((bc/ac)*ac-bc)> 10*%eps then pause,end
  159. //
  160. if abs(a*(a\b')-b')> 10*%eps then pause,end
  161. if abs(ac*(ac\b')-b')> 10*%eps then pause,end
  162. if abs(a*(a\bc')-bc')> 10*%eps then pause,end
  163. if abs(ac*(ac\bc')-bc')> 10*%eps then pause,end
  164. //
  165. [l u]=lu(a);
  166. if abs(l*u-a)> 10*%eps then pause,end
  167. [l u]=lu(ac);
  168. if abs(l*u-ac)> 10*%eps then pause,end
  169. //
  170. h1(5,5)=0;for k=1:5,for l=1:5, h1(k,l)=1/(k+l-1);end;end;
  171. if abs(inv(h1)-matrix('hilb',5))> 1.e-7 then pause,end
  172. //
  173. if abs(det(matrix('magic',5))-5070000)> 1.e-7 then pause,end
  174. //
  175. b=a*a';h=chol(b);
  176. if abs(h'*h-b)> 10*%eps then pause,end
  177. bc=triu(ac*ac');bc=bc+bc'-diag(real(diag(bc)));;h=chol(bc);
  178. if abs(h'*h-bc)> 10*%eps then pause,end
  179. //test de matqr
  180. a=rand(3,4);b=rand(3,4);ac=a+i*rand(3,4);bc=b+i*rand(3,4);
  181. //
  182. if abs(a*(1/a)*a-a)> 10*%eps then pause,end
  183. if abs(a*(i/a)*a-i*a)> 10*%eps then pause,end
  184. if abs(a*(a\ 1)*a-a)> 10*%eps then pause,end
  185. if abs(a*(a\ i)*a-i*a)> 10*%eps then pause,end
  186. if abs(ac*(1/ac)*ac-ac)> 10*%eps then pause,end
  187. if abs(ac*(i/ac)*ac-i*ac)> 10*%eps then pause,end
  188. if abs(ac*(ac\ 1)*ac-ac)> 10*%eps then pause,end
  189. if abs(ac*(ac\ i)*ac-i*ac)> 10*%eps then pause,end
  190. //
  191. if abs(a/b-a*(1/b))> 10*%eps then pause,end
  192. if abs(ac/b-ac*(1/b))> 10*%eps then pause,end
  193. if abs(a/bc-a*(1/bc))> 10*%eps then pause,end
  194. if abs(ac/bc-ac*(1/bc))> 10*%eps then pause,end
  195. //
  196. if abs(a\ b -(a\ 1)*b)> 10*%eps then pause,end
  197. if abs(a\ bc-(a\ 1)*bc)> 10*%eps then pause,end
  198. if abs(ac\ b-(ac\ 1)*b)> 10*%eps then pause,end
  199. if abs(ac\ bc-(ac\ 1)*bc)> 10*%eps then pause,end
  200. //
  201. //elemt-wise
  202. a=rand(3,2);ai=a+rand(3,2)*%i;
  203. de=2;
  204. if norm(ai.*de-ai*de,1) >1000*%eps then pause,end
  205. if norm(de.*ai-de*ai,1) >1000*%eps then pause,end
  206. if norm(ai.*de-ai*de,1) >1000*%eps then pause,end
  207. if norm(de.*ai-de*ai,1) >1000*%eps then pause,end
  208. de=de+3*%i;
  209. if norm(ai.*de-ai*de,1) >1000*%eps then pause,end
  210. if norm(de.*ai-de*ai,1) >1000*%eps then pause,end
  211. if norm(ai.*de-ai*de,1) >1000*%eps then pause,end
  212. if norm(de.*ai-de*ai,1) >1000*%eps then pause,end
  213. //
  214. de=2;def=de*ones(a);
  215. if norm(a./de-a./def,1) >1000*%eps then pause,end
  216. if norm(ai./de-ai./def,1) >1000*%eps then pause,end
  217. //
  218. de=2+3*%i;def=de*ones(a);
  219. if norm(a./de-a./def,1) >1000*%eps then pause,end
  220. if norm(ai./de-ai./def,1) >1000*%eps then pause,end
  221. de=2;def=de*ones(a);
  222. if norm(a.\de-a.\def,1) >1000*%eps then pause,end
  223. if norm(ai.\de-ai.\def,1) >1000*%eps then pause,end
  224. //
  225. de=2+3*%i;def=de*ones(a);
  226. if norm(a.\de-a.\def,1) >1000*%eps then pause,end
  227. if norm(ai.\de-ai.\def,1) >1000*%eps then pause,end
  228.  
  229. ///////////
  230. de=2;def=de*ones(a);
  231. if norm(de.\a-de.\a,1) >1000*%eps then pause,end
  232. if norm(de.\ai-def.\ai,1) >1000*%eps then pause,end
  233. //
  234. de=2+3*%i;def=de*ones(a);
  235. if norm(de.\a-def.\a,1) >1000*%eps then pause,end
  236. if norm(de.\ai-def.\ai,1) >1000*%eps then pause,end
  237. de=2;def=de*ones(a);
  238. if norm(de./a-def./a,1) >1000*%eps then pause,end
  239. if norm(de./ai-def./ai,1) >1000*%eps then pause,end
  240. //
  241. de=2+3*%i;def=de*ones(a);
  242. if norm(de./a-def./a,1) >1000*%eps then pause,end
  243. if norm(de./ai-def./ai,1) >1000*%eps then pause,end
  244. //
  245. [p,r]=qr(a);
  246. if abs(p*r-a)> 10*%eps then pause,end
  247. [p,r]=qr(a');
  248. if abs(p*r-a')> 10*%eps then pause,end
  249. [p,r,x]=qr(a);
  250. if abs(p*r*x'-a)> 10*%eps then pause,end
  251. [p,r]=qr(ac);
  252. if abs(p*r-ac)> 10*%eps then pause,end
  253. [p,r,x]=qr(ac);
  254. if abs(p*r-ac*x)> 10*%eps then pause,end
  255. //
  256. if abs(cond(diag([1 2 3 4]))-4)> 10*%eps then pause,end
  257. if abs(cond(diag([1 i 3 4]))-4)> 10*%eps then pause,end
  258. v=[1 2 3 4 5];
  259. if abs(norm(v,1)-15)> 10*%eps then pause,end
  260. if abs(norm(v,'inf')-5)> 10*%eps then pause,end
  261. if abs(norm(v,2)-sqrt(55))> 10*%eps then pause,end
  262. if abs(norm(v,'fro')-sqrt(55))> 10*%eps then pause,end
  263. v=v';
  264. if abs(norm(v,1)-15)> 10*%eps then pause,end
  265. if abs(norm(v,'inf')-5)> 10*%eps then pause,end
  266. if abs(norm(v,2)-sqrt(55))> 10*%eps then pause,end
  267. if abs(norm(v,'fro')-sqrt(55))> 10*%eps then pause,end
  268. v=[i 2 3 4 5];
  269. if abs(norm(v,'inf')-5)> 10*%eps then pause,end
  270. if abs(norm(v,2)-sqrt(55))> 10*%eps then pause,end
  271. if abs(norm(v,'fro')-sqrt(55))> 10*%eps then pause,end
  272. v=v';
  273. if abs(norm(v,1)-15)> 10*%eps then pause,end
  274. if abs(norm(v,'inf')-5)> 10*%eps then pause,end
  275. if abs(norm(v,2)-sqrt(55))> 10*%eps then pause,end
  276. if abs(norm(v,'fro')-sqrt(55))> 10*%eps then pause,end
  277. a=[diag([1 2 3]);[0 0 0]];
  278. if abs(norm(a,1)-3)> 10*%eps then pause,end
  279. if abs(norm(a,'inf')-3)> 10*%eps then pause,end
  280. if abs(norm(a,2)-3)> 10*%eps then pause,end
  281. if abs(norm(a,'fro')-sqrt(14))> 10*%eps then pause,end
  282. a=a';
  283. if abs(norm(a,1)-3)> 10*%eps then pause,end
  284. if abs(norm(a,'inf')-3)> 10*%eps then pause,end
  285. if abs(norm(a,2)-3)> 10*%eps then pause,end
  286. if abs(norm(a,'fro')-sqrt(14))> 10*%eps then pause,end
  287. a=[diag([i,2,3]),[0;0;0]];
  288. if abs(norm(a,'inf')-3)> 10*%eps then pause,end
  289. if abs(norm(a,2)-3)> 10*%eps then pause,end
  290. if abs(norm(a,'fro')-sqrt(14))> 10*%eps then pause,end
  291. a=a';
  292. if abs(norm(a,1)-3)> 10*%eps then pause,end
  293. if abs(norm(a,'inf')-3)> 10*%eps then pause,end
  294. if abs(norm(a,2)-3)> 10*%eps then pause,end
  295. if abs(norm(a,'fro')-sqrt(14))> 10*%eps then pause,end
  296. //
  297. a=rand(3,5);ac=a+i*rand(3,5);
  298. [u,s,v]=svd(a);u*s*v'-a;
  299. if abs(svd(a)-diag(s))> 10*%eps then pause,end
  300. [u,s,v]=svd(ac);u*s*v'-ac;
  301. if abs(svd(ac)-diag(s))> 10*%eps then pause,end
  302. //
  303. [u,s,v]=svd(a,0);u*s*v'-a;
  304. if abs(svd(a,0)-diag(s,0))> 10*%eps then pause,end
  305. [u,s,v]=svd(ac,0);u*s*v'-ac;
  306. if abs(svd(ac,0)-diag(s))> 10*%eps then pause,end
  307. a=a';ac=ac';
  308. a=rand(3,5);ac=a+i*rand(3,5);
  309. [u,s,v]=svd(a);u*s*v'-a;
  310. if abs(svd(a)-diag(s))> 10*%eps then pause,end
  311. [u,s,v]=svd(ac);u*s*v'-ac;
  312. if abs(svd(ac)-diag(s))> 10*%eps then pause,end
  313. //
  314. [u,s,v]=svd(a,0);u*s*v'-a;
  315. if abs(svd(a,0)-diag(s,0))> 10*%eps then pause,end
  316. [u,s,v]=svd(ac,0);u*s*v'-ac;
  317. if abs(svd(ac,0)-diag(s))> 10*%eps then pause,end
  318. //
  319. if abs(a*pinv(a)*a-a)> 10*%eps then pause,end
  320. if abs(ac*pinv(ac)*ac -ac)> 10*%eps then pause,end
  321. a=a';ac=ac';
  322. if abs(a*pinv(a)*a-a)> 10*%eps then pause,end
  323. if abs(ac*pinv(ac)*ac-ac)> 10*%eps then pause,end
  324. //
  325. if abs(rank(a)-3)> 10*%eps then pause,end
  326. if abs(rank(ac)-3)> 10*%eps then pause,end
  327. rand('seed',0)
  328. // matdsr  matdsc
  329. a=rand(4,4);ac=a+i*rand(4,4);t=a*a';tc=ac*ac';
  330. tc=triu(tc,1)+triu(tc,1)'+diag(real(diag(tc)));
  331. //
  332. //hess
  333. [u,h]=hess(a);
  334. if norm(u*h*u'-a,1) > 20*%eps then pause,end
  335. if abs(norm(h-hess(a),1))> 10*%eps then pause,end
  336. [u,h]=hess(ac);
  337. if norm(u*h*u'-ac,1) > 200*%eps then pause,end
  338. if abs(norm(h-hess(ac),1))> 10*%eps then pause,end
  339. [u,h]=hess(t);
  340. if norm(u*h*u'-t,1)  > 200*%eps then pause,end
  341. if abs(norm(h-hess(t),1))> 10*%eps then pause,end
  342. [u,h]=hess(tc);
  343. if norm(u*h*u'-tc,1)  > 200*%eps then pause,end
  344. if abs(norm(h-hess(tc),1))> 10*%eps then pause,end
  345. //schur
  346. [u,s]=schur(a);
  347. if norm(u*s*u'-a,1 ) > 200*%eps then pause,end
  348. if norm(s-schur(a),1 ) > 200*%eps then pause,end
  349. if norm(spec(a)-spec(s),1)> 200*%eps then pause,end
  350. [u,s]=schur(ac);
  351. if norm(u*s*u'-ac,1 ) > 200*%eps then pause,end
  352. if norm(s-schur(ac),1 ) > 200*%eps then pause,end
  353. if norm(spec(ac)-spec(s),1)> 200*%eps then pause,end
  354. [u,s]=schur(t);
  355. if norm(u*s*u'-t,1 ) > 200*%eps then pause,end
  356. if norm(s-schur(t),1 ) > 200*%eps then pause,end
  357. if norm(diag(s)-spec(t),1)> 200*%eps then pause,end
  358. [u,s]=schur(tc);
  359. if norm(u*s*u'-tc,1 ) > 200*%eps then pause,end
  360. if norm(s-schur(tc),1 ) > 200*%eps then pause,end
  361. if norm(diag(s)-spec(tc),1) > 200*%eps then pause,end
  362. // fonctions matricielles
  363. s=sqrt(t);
  364. if norm(t-s*s,1) > 200*%eps then pause,end
  365. s=log(t);
  366. if norm(t-exp(s)) > 200*%eps then pause,end
  367. s=sqrt(tc);
  368. if norm(tc-s*s,1) > 220*%eps then pause,end
  369. s=exp(tc);s=triu(s,1)+triu(s,1)'+diag(real(diag(s)));
  370. if norm(log(s)-tc,1)> 1.e-9 then pause,end
  371. if norm(sin(t)**2+cos(t)**2-eye,1) > 22*%eps then pause,end
  372. if norm(sin(tc)**2+cos(tc)**2-eye,1) > 20*%eps then pause,end
  373. //poly et root
  374. p=rand(5,1);pc=p+i*rand(5,1);x=poly(0,'x');
  375. if norm(sort(p )-sort(real(roots(poly(p,'x'))))) > 1000*%eps then pause,end
  376. if norm(sort(imag(pc))-sort(real(roots(poly(imag(pc),'x'))))) > 1000*%eps then pause,end
  377. //**
  378. if norm(t^(-1)-inv(t),1) > 200*%eps then pause,end
  379. if norm(tc**(-1)-inv(tc),1) > 200*%eps then pause,end 
  380. x=rand+i*rand;
  381. if norm((t**x)*(t^(-x))-eye,1) > 200*%eps then pause,end 
  382. if norm((tc**x)*(tc**(-x))-eye,1) > 200*%eps then pause,end
  383. //element op
  384. if norm(sin([0 pi/2 pi 3*pi/2])-[0 1 0 -1],1) > 10*%eps then pause,end
  385. if norm(cos([0 pi/2 pi 3*pi/2])-[1 0 -1 0],1) > 10*%eps then pause,end
  386. if norm(log([e 1/e e**2])-[1 -1 2],1) > 10*%eps then pause,end
  387. if norm(exp([1 -1 2])-[e 1/e e**2],1) > 10*%eps then pause,end
  388. if norm(sqrt([1 -1 4])-[1 i 2],1) > 10*%eps then pause,end
  389. if norm(sin(x)-(exp(i*x)-exp(-i*x))/(2*i),1) > 10*%eps then pause,end
  390. if norm(cos(x)-(exp(i*x)+exp(-i*x))/2,1) > 10*%eps then pause,end
  391. if norm(sqrt(x)^2-x,1) > 10*%eps then pause,end
  392. if norm(log(exp(x))-x,1) > 10*%eps then pause,end
  393. // triu
  394. z=poly(0,'z');a=[z 1 -z+1  8;z*z 10*z 5 -z;3 7 z+1 -1-z];
  395. [m,n]=size(a);mn=mini([m n]);a1=a;l=1;
  396. //
  397. for dg=-(m-1):0,
  398.    if norm(coeff(triu(a,dg)-a1),1)>%eps then dg,pause,end,
  399.    for k=1:l,a1(-dg+k,k)=0;end,l=mini([l+1,mn]);end
  400. for dg=1:n,
  401.    if norm(coeff(triu(a,dg)-a1),1)>%eps then dg,pause,end,
  402.    if dg>(n-m),l=l-1;end;for k=1:l,a1(k,dg+k)=0;end;end;
  403. //
  404. a=a';a1=a;[m,n]=size(a);mn=mini([m,n]);l=1;
  405. for dg=-(m-1):0,
  406.    if norm(coeff(triu(a,dg)-a1),1)>%eps then dg,pause,end,
  407.    for k=1:l,a1(-dg+k,k)=0;end,l=mini([l+1,mn]);end
  408. for dg=1:n,
  409.    if norm(coeff(triu(a,dg)-a1),1)>%eps then dg,pause,end,
  410.    if dg>(n-m),l=l-1;end;for k=1:l,a1(k,dg+k)=0;end;end;
  411. //
  412. a=a+%i*[1 2 3 4;5 6 7 8;9 10 11 12]';
  413. [m,n]=size(a);mn=mini([m n]);a1=a;l=1;
  414. //
  415. for dg=-(m-1):0,
  416.    if norm(coeff(triu(a,dg)-a1),1)>%eps then dg,pause,end,
  417.    for k=1:l,a1(-dg+k,k)=0;end,l=mini([l+1,mn]);end
  418. for dg=1:n,
  419.    if norm(coeff(triu(a,dg)-a1),1)>%eps then dg,pause,end,
  420.    if dg>(n-m),l=l-1;end;for k=1:l,a1(k,dg+k)=0;end;end;
  421. //
  422. a=a';a1=a;[m,n]=size(a);mn=mini([m,n]);l=1;
  423. for dg=-(m-1):0,
  424.    if norm(coeff(triu(a,dg)-a1),1)>%eps then dg,pause,end,
  425.    for k=1:l,a1(-dg+k,k)=0;end,l=mini([l+1,mn]);end
  426. for dg=1:n,
  427.    if norm(coeff(triu(a,dg)-a1),1)>%eps then dg,pause,end,
  428.    if dg>(n-m),l=l-1;end;for k=1:l,a1(k,dg+k)=0;end;end;
  429.  
  430. //test de clear
  431. x=[1,2,3];
  432. clear x,if exists('x')==1 then pause,end
  433. x=[1,2,3];y='afs';
  434. clear x y,if exists('x')==1| exists('y')==1 then pause,end
  435. X=[1,2,3];Y='afs';
  436. clear X Y,if exists('X')==1| exists('Y')==1 then pause,end
  437.